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The importance of triplet states in the photorelaxation dynamics of SO2 is studied by mixed quantum- 
classical dynamics simulations. Using the Surface Hopping including ARbitrary Couplings (Sharc) method, 
intersystem crossing processes caused by spin-orbit coupling are found occuring on an ultrafast time scale 
(below 100 fs) and thus competing with internal conversion. Comparing the dynamics including singlet and 
triplet states to the singlet-only dynamics, very different results are obtained for the populations of the 
respective states. However, the vibrational motion in the triplet manifold very much resembles the one in the 
singlet manifold. Consequently, the contribution of the triplet states may be difficult to detect in a broad 
range of experiments. 



I. INTRODUCTION 

Sulphur dioxide is an attractive molecular system 
which combines a very small system size with many dif- 
ferent excited states coupled in a complex way, offering 
a challenging photodynamics which can be studied with 
ultra-fast, time-resolved spectroscopic methods and mod- 
ern high-level ab initio techniques. 

In a companion paper,Sl (hereafter referred to as pa- 
per I) the experimental and theoretical efforts of the 
last 80 years to characterize the electronic excited states 
of SO9 have been surveyed in detail. Already in the 
1930s, ^ ^ high-resolution absorption spectra of SO2 were 
recorded. The absorption spectrum features three promi- 
nent bands in the near- and medium-UV range, named 
the forbidden band and the first and second allowed 
bands .El Especially the first allowed band between 3.6 eV 
and 5.1 eV was subject to numerous analyses (reviewed 
by Herzber^and Heicklen et al.l^ because of its intricate 
structure. From a modern point of view, this band sys- 
tem arises from the transition from the ground state to 
the vibronically coupled l^A2/l^Bi system (the two low- 
est excited singlet states). Miiller and Koppel^ treated 
this singlet state system using full-dimensional quantum 
dynamics (QD) and found that the system remains pri- 
marily on the lower adiabatic potential energy surface 
(PES) after excitation. An extension of this study has 
been published very recently.!^ 

However, there is significant evidenc^^H^I' that the ex- 
cited state dynamics within the first allowed band sys- 
tem is also affected by the presence of triplet states. SO2 
shows in the region between 3.1 eV and 3.6 eV a weak ab- 
sorption profile (the forbidden band), which was shown 
to arise from excitation to triplet states by Douglas^^ by 
means of the Zeeman effect. Even though the number, 
location and character of the triplet states could not be 
determined by early spectroscopic means, modern ab ini- 
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tio methods showed that three states of ^^2, ^Bi and 
symmetry are present at energies slightly lower than the 
corresponding singlet states. Consequently, intersystem 
crossing (ISC) between singlet and triplet states could be 
plausible and the photodynamics of SO2 might be influ- 
enced by both ISC and internal conversion (IC) processes. 

Very recently, both spectroscopicP^"^^ and 
theoretical^^ studies suggested that in a number 
of systems ISC indeed can compete with IC on an ultra- 
fast time scale. In order to provide more evidence for 
these new findings, several experimental methods based 
on high-resolution time-resolved spectroscopy have been 
developed. ^^ ^^ Further experimental indications for the 
case of SO2 can be found in paper I. 

In the present paper, we attempt to unravel the ex- 
cited state dynamics of SO2 theoretically. To this aim, 
ab initio surface hopping molecular dynamics (MD) is 
employed using the Sharc code.^^ Sharc (Surface Hop- 
ping including ARbitrary Couplings) can treat IC and 
ISC, mediated by non-adibatic couplings and spin-orbit 
couplings, respectively, on the same footing. This is 
achieved by performing surface hopping in a basis of 
spin-orbit-coupled electronic states, which are obtained 
by diagonalization of the potential energy matrix includ- 
ing the spin-orbit couplings. The simulations are fo- 
cused on the energy range corresponding to the forbid- 
den and first allowed band of the experimental absorption 
spectrum.!^ The applicability of Sharc on SO2 is vali- 
dated by running also trajectories on the singlet-manifold 
only, whose results are compared with those obtained by 
Miiller and Koppel using accurate quantum dynamical 
(QD) calculations.^ Other applications of Sharc can be 
found in Refs. 22. i3lH33l 



II. METHODOLOGY 

In the following, a brief summary of general surface 
hopping method as well as the specific Sharc variant is 
given. Further details on Sharc can be found in Ref. [29l 
A summary of the employed ab-initio methods and initial 
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conditions generation is also given below. 



A. Surface hopping 

The widely applied fewest switches criterion for surface 
hopping was devised by Tully in 1990.^"^ It allows incor- 
porating non-adiabatic effects into semi-classical ab initio 
dynamics, where the nuclei can be propagated on only 
one PES at a time, by means of state hops according to 
jumping probabilities. To this end, the electronic wave- 
function is expanded in the basis of adiabatic eigenstates 
IV^Q,) of the Hamiltonian: 



(1) 



The state coefRcients Cq, are obtained by numerical inte- 
gration of the time-dependent Schrodinger equation: 



d_ 
dt 
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■^c„[i(5^«^S+vT^«], (2) 



where v is the velocity vector, is the adiabatic energy 
of state a and 



(3) 



is the non-adiabatic coupling vector. 

Based on the state coefficients, the hopping probabil- 
ities from the current state into any other state can be 
computed according to: 
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K(c;c«[ii/^«+vT^«]). (4) 



B. Surface hopping including arbitrary couplings 




FIG. 1. Geometry of the molecule, internal coordinates and 
symmetry elements of SO2. The experimental values are 
ri=r2=1.432A and 6'=119.5°. 



so that the final equation for the evolution of the state 
coefficients is 



d_ 

m 



= - 



i(U^HU). 



(7) 



The motion of the nuclei is governed by gradients ob- 
tained numerically from the unperturbed gradients un- 
der the approximation of constant spin-orbit coupling el- 
ements, which is a good approximation as long as the 
coupling elements or their change are small. To date no 
analytical gradients for spin-orbit-perturbed states are 
available, but work is in progress. 

Using the described unitary transformation, the ob- 
tained states are not spin-pure states anymore because 
singlets and triplets mix via spin-orbit coupling. How- 
ever, in the basis of spin-orbit-coupled electronic states, 
IC and ISC are not fundamentally different anymore^ 
and can be conveniently investigated using non-adiabatic 
dynamics. 



In the scheme outlined above, the electronic Hamil- 
tonian was assumed to be diagonal in the basis of the 
wavefunctions IV^^). If perturbation terms are added 
to the Hamiltonian, off-diagonal elements between the 
states under consideration may be introduced: 



^(1) 



(0) 



Such perturbations include for example spin-orbit cou- 
pling, but also the interaction with laser fields. Within 
the Sharc algorithm, the Hamiltonian H^^^ is diagonal- 
ized by an unitary transformation. To obtain a consistent 
adiabatic picture, the non-adiabatic couplings are trans- 
formed into the new basis {IV^i^^)} according to: 



i) = (utT(°)u) +(UtVU^ 
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C. Ab initio methods 

Throughout this paper, the nuclear motion of SO2 is 
discussed in terms of a set of internal coordinates con- 
sisting of the average bond length, rgym, one half of the 



(5) bond length difference, Tasym, and the bond angle, 9: 
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where ri, r2 and are defined in figure [T] The MD 
calculations themselves were carried out in cartesian co- 
ordinates. 

The on-the-fiy electronic structure calculations re- 
quired for the MD simulations were performed with 
the MOLPRO 2010.1 package, which was interfaced 
to Sharc. In order to achieve a description of the 
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electronic states that is both feasible and accurate, we 
used the Complete Active Space Self-Consistent Field 
(CASSCF)P approach combined with the def2-SVFESl 
basis set to obtain preliminary wavefunctions of the 
three lowest-lying singlet and three lowest-lying triplet 
states. The employed active space comprises 12 elec- 
trons in 9 orbitals, so that this prescription is denoted 
as CASSCF(12,9)/def2-SVP. State-averaged calculations 
with equal weights were performed either on the three 
singlet states or the total six (three singlet and three 
triplet) electronic states. The so-obtained wavefunc- 
tions were subsequently correlated using Molpro's in- 
ternally contracted multireference configuration interac- 
tion (MRCI) formalism.EH The MRCI method was found 
to be necessary since the CASSCF potentials are other- 
wise too flat for large bond lengths - a fact already noted 
by Katagiri et al.^^ For the MRCI reference space, the 
three lowest active orbitals were constrained to contain 
at least 5 electrons, while the three highest active orbitals 
were constrained to a maximum of two electrons. To fur- 
ther speed up the calculations, configurations containing 
exactly one electron in the external orbital space were 
omitted. The spin-orbit matrix elements were calculated 
using the effective Fock-type spin-orbit operatoi^ imple- 
mented in MOLPRO. 



Gradient modifications 



each atom A were recomputed as follows: 
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In this way, the directions of the velocity vectors are not 
changed. The nuclear motion according to these gradi- 
ents was determined with the Velocity- Verlet algorithm^ 
using a time-step of 0.5 fs. 



E. Initial conditions 

The generation of the initial geometries and velocities 
was carried out with the Newton-X program suite . ^-^^ 
A quantum harmonic oscillator Wigner 

distributioiPEH 

was generated using numerical frequencies obtained with 
RI-MP2^^ and the basis set cc-pVTZ, as implemented in 
the TuRBOMOLE 6.2 suite. 

The obtained initial conditions displayed bond lengths 
between 1.36 A and 1.54 A and bond angles between 
110 ° and 128 °. The initial kinetic energies showed an 
average of 0.09 eV, with a maximum of 0.5 eV. 

All trajectories started in one of the excited states at 
t=0, which is equivalent to a delta pulse excitation. The 
initially populated state was determined stochastically 
based on computed transition probabilities Poa according 
to: 



Since Molpro is not able to calculate MRCI gradients 
analytically, approximate MRCI gradients were obtained 
from the CASSCF ones. The component of the CASSCF 
gradient in the direction of the trajectory's current mo- 
tion was substituted by a term obtained from the MRCI 
energy difference of the last and the current trajectory 
step. The new gradient reads as follows: 



g gCAS + 



AR AR 



AR AR 
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where gcAS is the gradient of the CASSCF energies, 
E^Ri) and £^(Ri_i) are the MRCI energies of the cur- 
rent and last step, respectively, and AR = R^ — R^-i 
is the displacement vector between the last and current 
geometry. 

In order to account for total energy conservation (de- 
spite small remaining inconsistencies), the kinetic energy 
is scaled after each step. To accomplish this, the kinetic 
energy of each atom A is recalculated according to: 
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Here, £^tot — ^pot is the kinetic energy required by en- 
ergy conversation. The kinetic energy is divided among 
the atoms proportionally to the unsealed kinetic energies. 
With the scaled kinetic energies, the velocity vectors of 



Poa — 



foa 

EL 



E m 

^0/3 



(12) 



where Eoa is the energy difference between ground state 
and state a and foa is the oscillator strength for this tran- 
sition. For the determination of the initial state, for each 
initial geometry one state in the energy band between 
3.5 eV and 5.5 eV was chosen based on the probabilities 
Poa- The state coefficient of state was set equal to 1 
and Cf3^a equal to 0. 



III. RESULTS 

A. State terminology 

To compare theoretical results to experimental find- 
ings, a note on nomenclature is in order, since spectro- 
scopic results are usually discussed in terms of diabatic 
states (where the wavefunction character is preserved), 
while ab initio MD simulations typically employ adia- 
batic states (where the potential coupling between states 
vanishes ).^^ 

Figure |2] shows the PESs of the lowest-lying five adi- 
abatic states that are relevant to the dynamics. Exem- 
plarily, the PES are shown in the (rasym,^) space. SO2 
exhibits C2V symmetry if Tasym = 0; otherwise, the point 
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FIG. 2. PESs of (a) singlet states 5*1 (lower surface) and & 
(upper surface) and (b) of triplet states Ti (lowest), T2 (inter- 
mediate), and T3 (highest). Colors indicate the wavefunction 
character of the adiabatic states, where purple corresponds to 
IM2 or 1^A2, cyan to l^Bi or l^Bi and green to 1^B2. For 
^asym 0, the dominating character is indicated by red (A2), 
blue (Bi) or grey (mixed). 

group changes Cg. Within C2V symmetry, all the con- 
sidered electronic states have different state symmetry 
and hence do not mix. The diabatic singlet states have 
state symmetry 1^742 and l^^i and are shown in purple 
and cyan, respectively, in figure [2] (a). For the remain- 
ing parts of the PESs, where rasym 7^ 0, both singlet 
states have symmetry A". Consequently, the two dia- 
batic states mix and give rise to the two adiabatic states 
Si (lower surface) and ^2 (upper surface) , which are used 
in the MD simulation. However, the contribution of the 
diabatic states to these mixtures depends on the geom- 
etry. We have indicated the predominant wavefunction 
character of the adiabatic (mixed) states at each point. 
Red corresponds to the character of the 1^742 state and 
blue to the one of l^^i. 

The same color scheme is used for the triplet states, 
shown in figure [2] (b) . Additionally, the state 1^52/1^^' 
is shown in green. Since even in Cs symmetry this state 
does not mix with the A" states, it retains its wavefunc- 
tion character along the whole surface. The adiabatic 
triplet states used in the MD simulations are called Ti, 
T2 and T3 and are obtained from l^A'', 2^A'' and l^A' 
by ordering the states strictly according to the energy at 
each geometry. For example, the T2 consists of the visi- 




FIG. 3. PESs of the two lowest excited singlet surfaces for 
^sym=l-5 A (a) and rsym=1.7 A (b). PESs of the three low- 
est triplet surfaces for rsym=l-5 A (c) and rsym=1.7 A (d). 
The ground state is given as a grey-scale contour map in 
panel (a) and the Franck-Condon point is indicated by an ar- 
row. Darker shades of blue (for singlets) or red (for triplets) 
indicate a stronger non-adiabatic coupling between the re- 
spective states. Singlet-triplet interaction regions (energy dif- 
ference <0.05 eV?^400 cm~^) are indicated in yellow. 



ble part of the green surface and the lower cone (pointing 
upwards) situated energetically above the green surface. 
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B. Potential energy surfaces 

In figure [sj panels (a) and (b) depict ttie potentials 
of the two excited singlet states in the same coordinate 
space as in figure |2] but using different axis ranges. In 
panel (a), rsym=l-5 A, which is similar to rgym at the 
ground state equilibrium geometry. The ground state 
energy is given as a greyscale contour to indicate the 
Franck-Condon region. In panel (b), rsym=l-7 A, which 
exemplarily represents the more stretched geometries 
found during the course of the MD simulation (see be- 
low). Regions of strong non-adiabatic coupling between 
the singlet states are shaded in dark blue, while yellow 
indicates regions where triplets are close to singlet states 
(less than 0.05 eV ^ 400 cm~^). In panels (c) and (d), 
the potentials of the three lowest-lying triplet states are 
shown analogously. 

In both panels (a) and (b), the singlet states are con- 
nected via a conical intersection (CI) directly situated 
above the Franck-Condon region (in (a) at ^=118.5 ° and 
rasym=0.0 A; in (b) at l9=120.0 ° and also rasym=0.0 A). 
The ^2 surface forms a narrow funnel without any in- 
tersection with triplet states. In contrast, the Si inter- 
sects with all three triplet states. Close to the CI, the Si 
crosses with the T3 , while Ti and T2 come close for larger 
values of Tasym- At large values for rgym, the Si surface 
exhibits a double- minimum (see panel (b)). 

Since the triplet states have basically the same wave- 
function character as their corresponding singlet states, 
the PESs are very similar. In panel (c), the T2/T3 CI is 
located at ^=111.0 °, rasym=0.0 A and in panel (d) at 
^=113.5 °, rasym=0.0 A. 



C. Spectrum 

Using the excitation energies Eoaj and the dipole mo- 
ments jj^oaj of the initial geometries of all trajectories j 
for all states a, an approximation to the absorption spec- 
trum has been simulated via a Gaussian convolution: 

g^{E) = X: foaj ■ exp (-^^^§^) • (13) 

Here, ga{E) is the part of the absorption spectrum aris- 
ing from the transition from the ground state to state a, 
depending on the energy E. The parameter e=0.027 eV 
describes the width of the Gaussian broadening. The 
spectrum is based on the properties of 280 initial condi- 
tions. 

The total absorption spectrum is the sum of ga {E) over 
all states a. To account for the first allowed band of SO2, 
only the two adiabatic states and ^2 need to be con- 
sidered. The obtained total spectrum together with its 
contributions from the transitions S'o ^ Si and ^ 5*2 , 
as well as the experimental spectrunPSl are depicted in 
figure [4j Both spectra were normalized. The agreement 




Energy (eV) 

FIG. 4. Absorption spectrum arising from excitation to 5*1 
and S2 (dark and light blue, respectively). The experimental 
spectrum of Golomb et al.^° is given in grey. 



between theory and experiment is satisfactory, consider- 
ing the overall energy range and the general shape of the 
band, although the calculated spectrum is slightly red- 
shifted. The shape of the high-energy fiank of the spec- 
trum is well reproduced, while the low-energy fiank gives 
too high oscillator strengths for energies below 3.9 eV. 
Naturally, the vibrational structure of the spectrum can- 
not be reproduced by this semi-classical method. 



D. Dynamics in the singlet-manifold 

In order to validate our method, we first performed 
a simulation of an ensemble of 42 trajectories, includ- 
ing only the three lowest-lying singlet states So to ^2. 
Thus, the results are directly comparable to those ob- 
tained by Miiller and Koppel, lousing full-dimensional QD 
on pre-calculated PESs for Si and 5*2. The PESs used in 
their simulations were obtained from the earlier work of 
Weis.l^^ 

Figure [5] shows the time-dependent populations of Si 
and ^2 as obtained by Miiller and Koppel (grey curves). 
At t=0, both states are populated, but absorbs a 
larger fraction of the population. During the first fem- 
toseconds, some small population fiuctuations take place. 
After approximately 10 fs, the 5*2 starts to transfer pop- 
ulation to the Si and gets completely depopulated after 
40 fs. After about 65 fs, again some population is trans- 
ferred from to ^2, because the wavepacket /ensemble 
returns to the interaction region of and ^2 near the 
Franck-Condon point. 

Figure [5] also depicts the time evolution of Si and 5*2 
obtained with Sharc. As it can be seen, the initial pop- 
ulation ratio obtained with Sharc is comparable to that 
predicted by QD. In the MD simulation, however, the 
population transfer from Si to S2 is faster (during the 
first 5 fs). Starting at about 15 fs, the 5*2 begins to con- 
tinuously transfer population to the Si until the former 



6 



Q. 
O 
CL 

> 
_rp 

CL> 



1.0 



0.8 



0.6 



0.4 



0.2 



0.0 




MP 

Si 

S2 



QD 



Si 

S2 



20 40 60 80 100 
Propagation Time (fs) 

FIG. 5. Relative population of 5*1 and & in the QD simulation 
of Miiller and Koppel (black and grey) and the corresponding 
results of an MD simulation including 42 trajectories (blue). 
Data taken from their paper. ^ 
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FIG. 6. Relative population of the six states in the dynamics 
simulation including triplet states. 



is completely depleted at around 60 fs. Between 70 and 
80 fs the 5*2 is repopulated, slightly delayed compared 
to the results of Miiller and Koppel. By the end of the 
propagations (100 fs), both methods predict a similar 
population exchange between both electronic states. In 
summary, despite different levels of theory in the QD and 
the MD simulations, a satisfactory agreement is found 
between both methods. 



E. Dynamics in the singlet-triplet manifold 

Encouraged by the general agreement found in the dy- 
namics performed within the singlet-manifold, an ensem- 
ble of 207 trajectories (171 starting in and 36 in S2) 
was propagated for 800 fs, now allowing the interaction 
with triplet states. In the following analysis, special em- 
phasis is put on the first 100 fs of the simulation. 

Figure [6] displays the relative populations of all singlet 
and triplet states depending on time in logarithmic scale. 
For simplicity, the populations of the three Ms compo- 
nents of each triplet state were summed up. The colors 
denoting the states in figure [6] (singlet states in shades of 
blue, triplets in shades of reel) are used consistently in all 
subsequent figures to allow for quick recognition of the 
respective state. 

In accordance with dynamics on the singlet-manifold, 
the initial population ratio is approximately 0.8:0.2 for 
Si:S2- The triplet states are unpopulated initially (note 
that figure |6] starts at t=l fs). However, the course of 
the dynamics shows remarkable differences with respect 
to the singlet-only case. Already after less than 10 fs, the 
population of Si and - to a lesser extent - ^2 decreases, 
while the one of T3 increases. The lower- lying Ti and 
T2 are populated only after around 30 fs. At the same 
time, the population of is diminished substantially. 
For later times, the states' populations fluctuate around 



some mean value: Si keeps about 50% of the total pop- 
ulation, Ti and T2 keep 20% each, and 5*2 and T3 the 
remainder. During the course of the simulation, also a 
small number of trajectories relaxed to the ground state 

In order to obtain a detailed interpretation of the state 
interactions, the numbers of trajectories hopping from 
state /3 to state a at time step t^, denoted AA^^Q,(t^), 
are analyzed. We use a Gaussian convolution in order to 
convert the discrete ANi^^itn) into more comprehensible, 
continuous data. The obtained "density of state hops" 
at time t is 



J^I3a{t) = AN(3a{tn) ' exp 



{tn - tf 

2r2 



(14) 



with a broadening parameter of r=5 fs. In the data so 
depicted, the area under the curves of M^a{t) is propor- 
tional to the number of state hops in a given time interval. 
The density of state hops is given in figure [7| for the first 
300 fs. 

Additionally, in order to extract time constants for 
the population transfer rate, we cumulate the number 
of state hops and fit all early transfer events separately 
using exponentials and all later (> 100 fs) according to 
zero-order kinetics (for an explanation of this distinction 
between early and late events, see below). 

Figure [7^ shows the interaction of 5*2 with Si. Di- 
rectly after the start of the simulation and again after 
approximately 35 fs, ^2 transfers population to ^i. For 
later times, ^2 does not show notable interaction with 
Si. The fitting procedure gave time constants of 7 and 
8 fs, respectively, for the two transfer events. Combined, 
16% of all trajectories were transferred to ^i. 

Ti gets populated from (panel (b)) after 30 fs and 
immediately afterwards the same amount of population 
is transferred back (boxes (2) and (3)). Consequently, 
during the first 100 fs, no net transfer from to Ti is 
observed. 
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FIG. 7. Density of net hops between interacting pairs of 
states. Dashed boxes enclose events discussed in the text. 
Colors denote the direction of the population transfer as de- 
noted in the labels. Common vertical scales have been used. 



After 20 fs, Si transfers a large number of trajectories 
to T2 (panel (c), box (4)) and again after 40 fs (box (5)). 
Treating the interval from 20 fs to 100 fs as one event, 
an effective time constant of 13 fs and a total transfer of 
33% of all trajectories is found. 

The interaction between states and T3 in panel (d) 
provides the earliest ISC process. The transfer of popu- 
lation to T3 (box (6)) starts immediately and during the 
first 100 fs, about 15% of all trajectories switch to T3 
with an effective time constant of 8 fs. 

The states Ti and T2 also interact with each other, 
see panel (e). Box (7) shows a transfer from T2 to Ti, 
starting after 20 fs. The exponential fit of the cumulated 
net hops yields a time constant of 11 fs. 

The last pair of states exhibiting notable interactions 
is T3 and T2, which are connected via a CI similar to the 
one linking S2 and Si. Panel (f) of figure [t] shows that 
the CI facilitates the ultrafast decay from T3 to T2 as 
soon as T3 is populated. This process is slightly slower 
than the ones already mentioned, with an estimated time 



constant of 17 fs. 

The state interactions for later times than 100 fs ex- 
hibit characteristics very different from the early part of 
the simulation. 5*2 and Ti interact significantly less with 
Si. However, states Si, T2 and T3 exchange population 
in a circular fashion, shown in the three boxes labeled 
(9). This cycle consists of three steps, T2 Si in panel 
(c), 5'i Ts in panel (d) and T3 T2 in panel (f). The 
cumulated net hops of these three interactions for times 
greater than 100 fs were fitted according to zero-order 
kinetics. Time constants of 10, 13 and again 13%/lOOfs 
(% of all trajectories per 100 fs) were found, respectively. 

Figure [8] shows the evolution of the internal coordi- 
nates of all trajectories in all states. In panel (a), the dy- 
namics in the asymmetric stretch mode is depicted. The 
lower-lying excited states Si , Ti and T2 all show a strong 
excitation of this mode. T2 and Ti get only populated 
once the trajectories has moved to large enough values 
of r. 



asym- 



All three states show almost exactly the same 
motion with regards to oscillation period and amplitude, 
thus we define these states as one group in the following 
discussion. On the contrary, a very different kind of mo- 
tion is exhibited by 5*2 and T3, which do not show any 
excitation of the asymmetric stretch mode. These two 
states therefore are classified as a second group. 

The classification of the excited states into the groups 
Si , Ti and T2 on the one hand and 5*2 and T3 on the other 
hand is confirmed by the other modes. For the symmetric 
stretch mode in panel (b), an oscillation period of 70 fs 
and maximum rgym of 1.8 A is found for the first group, 
while the second group performs one oscillation in 45 fs 
with a maximum rgym of slightly above 1.7 A. The bend- 
ing mode (panel (c)) of Si, Ti and T2 show a period of 
70 fs and a minimum angle of 95 °, the values of ^2 and 
Ts are 45 fs and 105 °. 



IV. DISCUSSION 

A. Potential energy surfaces 

The excitation in SO2 takes place from the l^Ai 
ground state to the bright l^Bi state. In the Franck- 
Condon region, this state crosses with 1^A2 and gives 
rise to a C I. Thus, in our adiabatic picture (see subsec- 
tion III A) the bright state l^^i contributes to both 
and S2 (for 6 < 118.5° mainly 6*2, otherwise mainly Si) 
and consequently both Si and ^2 are populated by the 
initial delta pulse. Because the CI is located right at the 
Franck-Condon region, IC takes place immediately after 
excitation. Another interesting feature of the PESs is 
the location of singlet-triplet interaction regions, favour- 
ing ISC (see figure [6|. Firstly, the intersection of Si and 
Ts circularly surrounds the S2/S1 CI and the Franck- 
Condon region. Therefore, all trajectories leaving the 
Franck-Condon region are forced to pass through this in- 
tersection. There, the spin-orbit coupling elements are 
on the order of 30 cm~^ and thus able to transfer a re- 
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FIG. 8. Evolution of the internal coordinates of all trajectories during the first 200 fs. 



markably large number of trajectories to the triplet PES. 
Secondly, the regions of proximity of S\ and both T\ 
and T2 coincide with the outer turning point of the 
asymmetric stretch mode, see Figure 3](b). Since at the 
turning point the potential is rather flat, the trajectories 
spend a long time in this region after excitation of the 
asymmetric stretch mode. Also in this region the spin- 
orbit couplings amount to about 30 cm~^ and allow for a 
noteworthy ISC yield. Since in this part of the PES 
Ti and T2 all become very close in energy, simultaneous 
ISC to Ti and T2 becomes possible. 



B. Spectrum 

As can be seen in Figure |4j both transitions 5*0 Si 
and ^0 S2 contribute to the spectrum to a similar 
extent. The resulting calculated absorption spectrum is 
in reasonable agreement with the experimental one, al- 
though small deviations are patent. First, the calculated 
overall band shape is slightly red-shifted, by approxi- 
mately 0.1 eV. This may be attributed to the double-zeta 
quality of the basis set since preliminary calculations in- 
dicate that larger basis sets give slightly higher excita- 
tion energies for all considered excited states. Despite its 
lesser accuracy, the double-zeta basis set was employed 
for the sake of computational performance during the dy- 
namics. 

Second, while the high-energy slopes of both spectra 
are quite similar, the low-energy slope of the calculated 
spectrum features too high intensities. In the experimen- 
tal spectrum, at around 3.8 eV the absorption is already 
close to zero, while in the calculated spectrum the slope 



extends until 3.5 eV. It is noteworthy that the experi- 
mental 0-0 transjiion for this band has been located at 
around 3.46 eV,*^so the obtained excitation energies are 
plausible and only the intensities are overestimated. 

Finally, it is obvious that the vibrational structure of 
the spectrum is not reproduced by our simulation. A de- 
scription of this phenomenon goes beyond of the scope of 
this work, since it requires a quantum-mechanical treat- 
ment of the nuclear motion. A semi-classical simulation 
cannot account for quantization of the vibrational energy 
levels and also is not able to deliver nuclear wavefunction 
overlaps, which are necessary to obtain Franck-Condon 
factors. Thus, a correct calculation of the total absorp- 
tion is not possible within this approach. The absence 
of the Franck-Condon factors might also explain the too 
intense low-energy slope in the calculated spectrum. 

In conclusion, given the employed theoretical model, 
the simulated spectrum is in good agreement with the 
experiment and confirms the correctness of the initial 
conditions. 



C. Dynamics in the singlet-manifold 

A good agreement is obtained between the QD simu- 
lation by Miiller and Koppel^ and our singlet-only MD 
simulations (recall Figure [5|. Right after excitation, the 
states and 5*2 briefly interact while still in proximity 
to the CI. However, already after ca 10 fs the popula- 
tion on Si has moved sufficiently far away from the CI 
to cease Si S2 population transfer. Afterwards, only 
the 5*2 population stays close to the interaction region, 
since the cone focusses the population on the intersec- 
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tion. Consequently, between 10 fs and 60 fs a decay from 
5*2 to Si is found. At 60 fs, the 5*2 surface is completely 
depopulated in both QD and MD simulations. 

The use of C ASSCF in the QD simulation and of MRCI 
in the MD simulation lead to differences in the potentials' 
slopes, non-adiabatic coupling strengths and the exact 
location of the CI. Already Katagiri et al.^^ noticed that 
the CASSCF method gives unsatisfactorily flat potentials 
of the Si state near the O— S---0 dissociation limit. On 
the contrary, the dissociation energy as obtained on the 
MRCI level of theory (about 5.8 eV) is very close to the 
experimental one at 5.65 eV. The small discrepancies in 
the populations of the two simulations can be explained 
in this way. 

After 70 fs (QD) or 80 fs (MD), the population on the 
Si surface returns to the S1/S2 intersection area in the 
Franck-Condon region. At this point, the ^2 is partly 
repopulated, though the initial population is not reached 
in both simulation types. This is reasonable in the sense 
that the ensemble already exhibits a considerable spread 
and the CI is missed by part of the ensemble. 

In both studies, the wavepacket /ensemble basically 
moves all the time in Si^ despite having sufficient en- 
ergy to reach the 5*2 . This behaviour can be attributed 
to the wavepacket /ensemble avoiding the CI because 
of the surrounding double-minimum potential (see fig- 
ure [3|. Following the potential energy gradient, the 
wavepacket /ensemble moves into regions of ri 7^ r2, 
where the two surfaces avoid each other and the non- 
adiabatic coupling is small. Also, both studies show that 
for times after 80 fs a coherent motion is virtually lost. 
Note that the experiment (see paper I) indicates a long- 
lived coherent motion. Although we could only specu- 
late about the reasons for this disagreement in the QD 
case, we attribute the discrepancies between experiment 
and our MD results to our approximate MRCI gradients 
(see section II D[ ). We are aware of the shortcomings of 
these gradients, however, MRCI potentials are necessary 
to qualitatively the SO2 dynamics and these gradients 
are the best approximation currently available compati- 
ble with the Sharc methodology. 

To sum up the findings so far, it is shown that the 
semi-classical dynamics gives results in satisfying agree- 
ment with the quantum dynamical results of Miiller and 
Koppel. Both show a high population of Si at all times, 
with a fraction of the population oscillating between the 
two surfaces. 



D. Dynamics in the singlet-triplet manifold 

One of the most important results of this study is that 
the inclusion of triplet states leads to a completely dif- 
ferent picture than that obtained in the singlet-only dy- 
namics. ISC plays a significant role in the deactivation of 
SO2 and, as can be seen from figure [6| ISC even competes 
with IC on a timescale of tens of femtoseconds. Exper- 
imental results of paper I (figure 6 (d)) also hint at the 



participation of ISC in the ultrafast dynamics of SO2. 

The differences between the singlet-only and the 
singlet-triplet dynamics are discussed in detail in the fol- 
lowing. Already during the first 5 fs, in the singlet-only 
dynamics there is considerable population transfer be- 
tween Si and ^2 close to the CI, while in the singlet- 
triplet simulation this exchange is almost absent. In- 
stead, the Si population approaches the intersection with 
T3 very fast and a large number of trajectories hop to the 
triplet surface. 

The presence of T3 also governs the fate of ^2 at later 
times. In the singlet-only simulation, S2 is repopulated 
each time Si returns to the Franck-Condon region since 
the S2/S1 CI located there. In the singlet-triplet simula- 
tion, the Si/Ts intersection completely surrounds the CI 
and all trajectories returning to the Franck-Condon re- 
gion have to pass through the singlet-triplet intersection. 
This leads to T3 scavenging a large amount of incom- 
ing trajectories which could otherwise repopulate S2' It 
is mainly this mechanism which prevents S2 from con- 
tributing significantly at later times. 

However, the inclusion of the triplet states has its 
greatest impact on the dynamics on Si. This state acts 
as the central hub in the dynamics and directly interacts 
with all other excited states under consideration. There- 
fore, in contrast to the singlet-only case, the Si is heavily 
depleted by ISC. The large number of trajectories show- 
ing ISC can be rationalized by the fact that the Si in- 
teracts with the triplets over a wide range of geometries, 
while the singlet-singlet interaction is significant only in 
a very localized area (the CI). 

The analysis of the net hops between the surfaces, pre- 
sented in Figure [7| has been summarized in a global 
scheme (Figure |9| that collects the processes encountered 
in the photorelaxation dynamics of SO2. Each state is 
represented by a box, indicating also its respective initial 
population. The So is absent from the figure, since re- 
laxation to the ground state is unlikely in the time scales 
we are considering here. 

The population transfer events have been divided into 
two sets. In the first set - the early transfers happen- 
ing during the first 100 fs - mostly temporally localized 
events occur (solid arrows in the figure), while for later 
times we see the continuous population transfers (dashed 
arrows) of the second set. The temporal localization is 
caused by the trajectories moving collectively and pass- 
ing jointly through the different interaction regions. For 
each localized transfer event, the effective time constant 
r, the onset of this process and the percentage of the to- 
tal population transferred in this event is given next to 
the corresponding solid arrow. Only major events are in- 
dicated. It is immediately obvious that Si interacts with 
all other excited states. It can also be seen that ISC and 
IC start instantly and consequently, both mechanisms act 
on an ultrafast timescale. 

After approximately 100 fs, the ensemble of trajecto- 
ries is considerably spread across the surfaces and there- 
fore temporally localized events are not observed any- 
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FIG. 9. Overview over transfer processes. The state boxes 
contain the initial relative population. Solid arrows denote 
the population transfers within the first 100 fs. The labels of 
these arrows give the rate constant r, the onset (in bold) and 
the total population transferred during this process. Dashed 
arrows denote population transfer after the first 100 fs, where 
the labels give the rate of transfer. 



more. Instead, we find a continuous transfer of popula- 
tion among a subset of the considered states. Circular 
population transfer between Si, T3 and T2 leads to a 
dynamical equilibrium. To a lesser extent, Ti also con- 
tributes to the dynamical equilibrium by the alternative 
cycle S'l ^ T3 ^ T2 ^ Ti ^ 

In view of the occurrence of localized transfer events 
and the mentioned circular transfer, a description of the 
whole dynamics by means of overall exponential fits is 
not suitable in the case of the relaxation dynamics in 

SO2. 

Another very interesting feature of SO2, already noted 
in subsection |III E| with regard to figure [8) is the division 
of the considered PESs into two classes of PESs. States 
^i, Ti and T2 are characterized by a strong excitation 
of the asymmetric stretch mode because of their double- 
minimum potentials. On the contrary, the PESs of S2 
and T3 do not allow any asymmetric stretch. What is 
most interesting about these two observations is that the 
nuclear motion, i.e. the vibrational frequency, does not 
allow for the identification of the populated electronic 
state. This fact may help to explain the experimen- 
tal difficult}/^ ^ to observe the triplet states throughout 
decades of spectroscopic research. 



In states Si, Ti and T2, some trajectories show asym- 
metric stretch oscillations (figure [s] (a)) exclusively in 
one of the wells of the double-minimum potential (from 
^asym=0.0 to arouud 0.4 A and back) with a period of 
about 70 fs. Another portion of the ensemble overcomes 
the central barrier and oscillates in both wells (from 
^asym=-0-4 to -|-0.4 A), with a full-cycle period of 140 fs. 
Interestingly, the oscillation period of 70 fs is also found 
in the symmetric stretch and bending modes (figure [s] 
(b) and (c)). Accordingly, the ensemble returns to the 
Franck- Condon region after this time interval and simi- 
larly at multiples of 70 fs. 

From the spac ing o f the Clement's bands in the ab- 
sorption spectrunP^ (approximately 28 meV) one would 
expect a recurrence time of 148 fs. This time coincides 
with the asymmetric oscillation period across both po- 
tential wells found in the simulation, but is about double 
the first recurrence time of the ensemble in the Franck- 
Condon region. Our recurrence time of about 70 fs would 
lead to a spacing of 59 meV, which could in principle 
be hidden under the peaks spaced by 28 meV. In order 
to recover the relative experimental peak heights in the 
spectrum, the second recurrence needs to be more pro- 
nounced than the first one. We do not observe such a 
behavior. The reason is possibly that the ensemble is 
already too strongly dispersed prior to the second recur- 
rence due to our approximate gradients. We believe that 
the oscillation period of 70 fs is correct, but is only due to 
a partial recurrence to the Franck- Condon region. Con- 
sequently, the expected peaks in the TRPEPICO (time- 
resolved photoelectron photoion coincidence) spectra of 
paper I might be hidden under the peaks correspond- 
ing to the 140 fs period. The full recurrence after these 
140 fs is underestimated in our dynamics most proba- 
bly because of the approximate gradients. Our rationale 
is supported by the findings of Leveque et alJ published 
very recently. In their remarkable study, they discuss the 
role of the different spectroscopic states in detail and are 
able to achieve almost quantitative agreement with the 
experimental spectrum for the first time. 



V. CONCLUSION 

This work is - to the best of our knowledge - the first 
ab initio-based dynamics study to investigate the singlet- 
triplet interactions in sulphur dioxide. The comparison 
with singlet-only results from QD and ab initio MD re- 
veals that the inclusion of the triplet states is essential 
in order to understand the photodynamics of this sys- 
tem. Due to the presence of singlet-triplet intersections 
at accessible sites on the PESs, not only ultrafast ISC 
processes have been found, but those strongly influence 
the singlet-singlet interaction. 

After the initial delta pulse excitation, 80% of all tra- 
jectories start in the Si state and 20% in 6*2 . While 5*2 
quickly decays to Si via a S2/S1 CI, in Si the asym- 
metric stretch mode is strongly excited. This stretch- 
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ing of the S— O bonds leads the trajectories inevitably 
through the intersection regions Si/Ts and slightly later 
S1/T1/T2. Spin-orbit couplings on the order of 30 cm~^ 
are responsible for an efficient and ultrafast ISC in the re- 
gions of near-degeneracy. ISC together with IC between 
the triplet states is responsible for the establishment of a 
dynamical equilibrium after 100 fs. There, Si, T3 and T2 
exchange population circularly, while ^2 and Ti remain 
mostly static. The occurrence of ISC in our MD simula- 
tion fits nicely with the observation of signals attributed 
to quartett states in decay associated spectra in paper I. 

It was also found that the PESs of Si, Ti and T2 on 
one hand and 5*2 and T3 on the other hand are suffi- 
ciently similar to cause indistinguishable molecular mo- 
tion on these surfaces. This may explain the difficulties 
encountered in spectroscopic experiments to identify the 
participation of triplet states in the excited-state dynam- 
ics of SO2. We observed a dominant oscillation period of 
70 fs in our simulation, which was related to the spacing 
of the Clement's bands in the experimental absorption 
spectrum and the TRPEPICO spectra in paper I. 

Thus we have shown that our recently developed ab 
initio molecular dynamics program SharcP^ can provide 
significant insight into complex molecular processes in- 
cluding states of different multiplicities. 
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